/*
 *  SMImage.c
 *  Seams
 *
 *  Created by Charles Osmer on 4/10/09.
 *  Copyright 2009 __MyCompanyName__. All rights reserved.
 *
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <assert.h>
#include "SMImage.h"

#pragma mark Image Lifespan

const SMImage SMImage_Initializer = {NULL, 0, 0, 0};

// @ Create a new SMImage object.
// @    Supports 1, 2, or 3 planes
// @    All pixels are set to zero
void SMImage_Alloc(SMImage *image, const int nPlanes, const int width, const int height)
{
    assert(image != NULL);
    assert(nPlanes >= 1 && nPlanes <= 3 && width > 0 && height > 0);

    image->nPlanes = nPlanes;
    image->width = width;
    image->height = height;

    image->planes[0] = image->planes[1] = image->planes[2] = NULL;

    for (int p=0; p < nPlanes; p++) {
        image->planes[p] = malloc(height*sizeof(SMPlaneType *));

        for (int y=0; y < height; y++)
            image->planes[p][y] = calloc(width, sizeof(SMPlaneType));
    }
}

// @ Free all memory allocated to SMImage object.
void SMImage_Free(SMImage *image)
{
    SMImageAssert(image);

    for (int p=0; p < image->nPlanes; p++) {
        SMImagePlane plane = image->planes[p];

        for (int y=0; y < image->height; y++)
            free(plane[y]);

        free(plane);
        plane = NULL;
    }

    image->nPlanes = image->width = image->height = 0;
}

void SMImage_Copy(SMImage *copy, const SMImage *original)
{
    SMImageAssert(original);
    assert(copy != NULL);

    SMImage_Alloc(copy, original->nPlanes, original->width, original->height);

    for (int p=0; p < original->nPlanes; p++) {        
        for (int y=0; y < original->height; y++)
            memcpy(copy->planes[p][y], original->planes[p][y], original->width*sizeof(SMPlaneType));
    }
}

#pragma mark -
#pragma mark Colorspace

// @ Creates a new gray scale image from an RGB image
// @    Averages RGB values
SMImage SMImage_RGB2Gray(const SMImage *image)
{
    SMImageAssert(image);
    assert(image->nPlanes == 3);

    SMImage gray = SMImage_Initializer;
    SMImage_Alloc(&gray, 1, image->width, image->height);

    for (int y=0; y < image->height; y++) {
        for (int x=0; x < image->width; x++) {
            const int pixel = (image->planes[R][y][x]+image->planes[G][y][x]+image->planes[B][y][x]);
            gray.planes[GRAY][y][x] = pixel/3;
        }
    }

    return gray;
}

#pragma mark -
#pragma mark Point Operations

// @ Normalize histogram to [0 high]
void SMImage_Normalize(const SMImage *image, unsigned int high)
{
    SMImageAssert(image);
    assert(high > 0);

    const int w = image->width;
    const int h = image->height;

    int min = image->planes[0][0][0];
    int max = image->planes[0][0][0];

    for (int p=0; p < image->nPlanes; p++) {
        for (int y=0; y < h; y++) {
            for (int x=0; x < w; x++) {
                if (image->planes[p][y][x] < min) min = image->planes[p][y][x];
                if (image->planes[p][y][x] > max) max = image->planes[p][y][x];
            }
        }
    }

    for (int p=0; p < image->nPlanes; p++) {
        for (int y=0; y < h; y++) {
            for (int x=0; x < w; x++) {
                image->planes[p][y][x] = high*(image->planes[p][y][x]-min)/(max-min);
            }
        }
    }
}

#pragma mark -
#pragma mark Edge Dectection

SMImage SMImage_SobelEdges(const SMImage *image)
{
    SMImageAssert(image);

    // Sobel operators
	const int hx[3][3] = { 
        { -1, 0, 1 },
        { -2, 0, 2 },
        { -1, 0, 1 }};
    
	const int hy[3][3] = {
        {  1,  2,  1 },
        {  0,  0,  0 },
        { -1, -2, -1 }};

    SMImage edges = SMImage_Initializer;
    SMImage_Alloc(&edges, image->nPlanes, image->width, image->height);

    for (int p=0; p < image->nPlanes; p++) {        
        for (int y=1; y < image->height-1; y++) {
            for (int x=1; x < image->width-1; x++) {
                int gx = 0, gy = 0;

                for (int dy=0; dy < 3; dy++) {
                    for (int dx=0; dx < 3; dx++) {
                        const int pixel  = image->planes[p][y+dy-1][x+dx-1];
                        
                        gx += pixel*hx[dy][dx];
                        gy += pixel*hy[dy][dx];
                    }
                }
                
                edges.planes[p][y][x] = abs(gx) + abs(gy);
            }
        }
    }

    return edges;
}

#pragma mark -
#pragma mark Differential Synthesis
#if 0
SMImage SMImage_Laplacian(const SMImage *image)
{
    SMImageAssert(image);

    SMImage laplacian = SMImage_Initializer;
    SMImage_Alloc(&laplacian, image->nPlanes, image->width-2, image->height-2);

    // Laplacian operator
	const int h[3][3] = { 
        {  0, -1,  0 },
        { -1,  4, -1 },
        {  0, -1,  0 }};

    for (int p=0; p < image->nPlanes; p++) {
        for (int y=1; y < image->height-1; y++) {
            for (int x=1; x < image->width-1; x++) {
                int g = 0;

                for (int dy=0; dy < 3; dy++) {
                    for (int dx=0; dx < 3; dx++) {
                        const int pixel = image->planes[p][y+dy-1][x+dx-1];

                        g += pixel*h[dy][dx];
                    }
                }

                laplacian.planes[p][y-1][x-1] = g;
            }
        }
    }

    return laplacian;
}

SMImage SMImage_PoissionReconstruction(const SMImage *image, const SMImage *laplacian)
{
    SMImageAssert(image);
    SMImageAssert(laplacian);
    assert(image->nPlanes == laplacian->nPlanes);

    SMImage poisson = SMImage_Initializer;
    SMImage_Alloc(&poisson, image->nPlanes, laplacian->width, laplacian->height);
    
    for (int p=0; p < image->nPlanes; p++) {
        SMImagePlane plane = poisson.planes[p];

        for (int y=0; y < poisson.height; y++) {
            for (int x=0; x < poisson.width; x++) {
                plane[y][x] = image->planes[p][y][x];
            }
        }

        for (int i=0; i < 1000; i++) {
            for (int y=0; y < laplacian->height; y++) {
                for (int x=0; x < laplacian->width; x++) {       
                    int f = laplacian->planes[p][y][x];

                    int n = 0;
                    int m = 0;

                    if (y != laplacian->height-1) {
                        n += plane[y+1][x];
                        ++m;
                    }

                    if (x != laplacian->width-1) {
                        n += plane[y][x+1];
                        ++m;
                    }

                    if (y != 0) {
                        n += plane[y-1][x];
                        ++m;
                    }

                    if (x != 0) {
                        n += plane[y][x-1];
                        ++m;
                    }

                    plane[y][x] = (n+f)/m;
                }
            }
        }
    }

    return poisson;
}
#endif
